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Abstract 

The pancreatic islets of Langerhans are multicellular micro-organs integral to maintaining glucose homeostasis through 
secretion of the hormone insulin, p-cells within the islet exist as a highly coupled electrical network which coordinates 
electrical activity and insulin release at high glucose, but leads to global suppression at basal glucose. Despite its 
importance, how network dynamics generate this emergent binary on/off behavior remains to be elucidated. Previous work 
has suggested that a small threshold of quiescent cells is able to suppress the entire network. By modeling the islet as a 
Boolean network, we predicted a phase-transition between globally active and inactive states would emerge near this 
threshold number of cells, indicative of critical behavior. This was tested using islets with an inducible-expression mutation 
which renders defined numbers of cells electrically inactive, together with pharmacological modulation of electrical activity. 
This was combined with real-time imaging of intracellular free-calcium activity [Ca^T, and measurement of physiological 
parameters in mice. As the number of inexcitable cells was increased beyond ~15%, a phase-transition in islet activity 
occurred, switching from globally active wild-type behavior to global quiescence. This phase-transition was also seen in 
insulin secretion and blood glucose, indicating physiological impact. This behavior was reproduced in a multicellular 
dynamical model suggesting critical behavior in the islet may obey general properties of coupled heterogeneous networks. 
This study represents the first detailed explanation for how the islet facilitates inhibitory activity in spite of a heterogeneous 
cell population, as well as the role this plays in diabetes and its reversal. We further explain how islets utilize this critical 
behavior to leverage cellular heterogeneity and coordinate a robust insulin response with high dynamic range. These 
findings also give new insight into emergent multicellular dynamics in general which are applicable to many coupled 
physiological systems, specifically where inhibitory dynamics result from coupled networks. 
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introduction 

Most biological systems exist as dynamic multicellular structures 
where distinct functionalities are generated through cellular 
interactions. While important for proper function, the complexity 
in network architecture, cellular dynamics, as well as the presence 
of heterogeneity, noise and biological variability make the overall 
function of multicellular structures difficult to understand. 
Approaches to understanding coupled dynamical systems have 
handled this complexity by explaining system structure and 
function individually [1,2]. These two aspects are both of central 
importance when it comes to understanding the way living systems 
are organized and how their anatomy supports their function. 
Therefore, by employing network theory to inform or predict the 
architectural aspects of dynamical system models, we can better 
understand how structural properties can impact functional 
behaviors. One living system exhibiting complex multicellular 
dynamics, yet with a scale tractable for study with these 



approaches, is the islet of Langerhans where dysfunction generally 
leads to diabetes. As such the islet provides a physiologically 
relevant system in which we can examine properties of multicel- 
lular dynamical systems and discover behavior that is broadly 
apphcable. 

The islets of Langerhans are multicellular micro-organs located 
in the pancreas which maintain glucose homeostasis through the 
secretion of hormones such as insulin. Glucose-stimulated insulin 
secretion (GSIS) from P-cells within the islet is driven by glucose- 
dependent electrical activity. The metabohsm of glucose and 
increased ATP/ADP ratio inhibits ATP-sensitive K""" (Katp) 
channels, causing membrane depolarization. Activation of volt- 
age-dependent Ca"''" channels elevates intracellular free-calcium 
activity ([Ca^"^;) to trigger insulin granule exocytosis [3,4]. Defects 
at several points in this signaling pathway, including the Katp 
channel, can cause or enhance the risk of developing diabetes [5— 
8]. Despite the importance of this pathway, it is important to 
recognize P-cells do not act autonomously. Rather, like many 
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Author Summary 

As science has successfully broken down the elennents of 
many biological systems, the network dynamics of large- 
scale cellular interactions has emerged as a new frontier. 
One way to understand how dynamical elements within 
large networks behave collectively is via mathematical 
modeling. Diabetes, which is of increasing international 
concern, is commonly caused by a deterioration of these 
complex dynamics in a highly coupled micro-organ called 
the islet of Langerhans. Therefore, if we are to understand 
diabetes and how to treat it, we must understand how 
coupling affects ensemble dynamics. While the role of 
network connectivity in islet excitation under stimulatory 
conditions has been well studied, how connectivity also 
suppresses activity under fasting conditions remains to be 
elucidated. Here we use two network models of islet 
connectivity to investigate this process. Using genetically 
altered islets and pharmacological treatments, we show 
how suppression of islet activity is solely dependent on a 
threshold number of inactive cells. We found that the islet 
exhibits critical behavior in the threshold region, rapidly 
transitioning from global activity to inactivity. We there- 
fore propose how the islet and multicellular systems in 
general can generate a robust stimulated response from a 
heterogeneous cell population. 



tissues, there are extensive cell-cell interactions within the islet that 
govern overall function. For example, isolated P-cells exhibit 
heterogeneous sensitivities to glucose with a low overall dynamic 
range of GSIS [9-11], yet P-cells within the islet robustly release 
insulin. Connexin36 (Cx36) gap junctions mediate the electrical 
coupling between P-cells [12-14] which coordinates oscillations in 
electrical activity and insulin release across the islet, enhancing the 
pulsatile release of insulin and glucose homeostasis [13-15]. In the 
absence of coupling many cells in the islet also show spontaneous 
elevations in [Ca"^ ;; likely as a result of heterogeneities in glucose 
sensitivity [10,16]. Therefore, another equally important role gap 
junctions play is to coordinate a suppression of spontaneous 
electrical activity at lower glucose levels [17]. Given that basal 
regulation is integral to glucose homeostasis, electrical coupling 
and the coordinated electrical dynamics are a critical factor in the 
regulation of islet function and in diabetes. 

Multicellular electrical dynamics in the islet have been described 
as functional networks where synchronized changes in [Ca^"*]; 
indicate functional connectivity between cells [14,18,19]. Such 
network analysis has been applied to examine the dependence of 
[Ca^"*]; dynamics on the level of coupling and its regulation, and 
has indicated that P-cell connectivity is non-homogeneous with a 
small subset of connections dominating synchronized behavior. As 
part of this analysis, the network of functional connectivity can be 
approximated by a Boolean network which quantitatively 
describes changes in multicellular behavior, including changes in 
coupling strength, network size or network shape [14,20-22]. 
These studies have generally focused on the synchronization of 
[Ca^"*]; oscillations, and such synchronized oscillatory/pulsatile 
behavior has been similarly examined in other physiological 
multicellular systems [23-25]. However, few studies have theoret- 
ically examined the suppressive effect of electrical coupling in the 
islet and its ability to shape the glucose-regulation of electrical 
activity. This is particularly warranted given a recent study that 
showed how severe diabetes caused by expression of mutant Katp 
channels could be prevented through a modulation in gap junction 
coupling [26]. Therefore, details for how the network structure 



and composition facilitate a highly sensitive and robust response 
from a heterogeneous cell population remain to be determined. 

In this study we examine how electrical coupling within P-ceU 
networks in the islet provide resilience against heterogeneous cell 
populations to generate robust network responses. We first develop 
quantitative predictions derived from a Boolean approximation of 
the P-ceU network, where the dependence of [Ca^"^; on variations 
in the constituent cellular excitability and coupling is described. 
We then test these predictions using two experimental systems 
involving transgenic mice that express mutant Katp channels with 
increased or decreased ATP-sensitivity [27,28]. This creates 
defined populations of cells within the islet which are 'excitable' 
or 'inexcitable', and can be further used to examine how our 
theoretical predictions and experimental data extend to physio- 
logical regulation of glucose homeostasis. We next link the static 
Boolean network model predictions and experimental findings 
with a dynamic multicellular model of the islet which incorporates 
recent understanding of P-ceU electrophysiology [29,30]. We 
finally extend these experimental and theoretical measurements to 
a general case with a continuum of heterogeneous cellular 
behavior. 

A consistent feature in this study is the emergence of critical 
behavior as a result of P-ceU electrical coupling, where the islet 
exhibits a phase transition between globally active and inactive 
states as cellular excitability approaches a critical threshold value. 
We discuss how the robust functionality that emerges at the 
multicellular level is not only relevant to the islet of Langerhans 
and its dysfunction in diabetes, but also to the function of other 
multicellular biological systems. 

Results 

Boolean Network Model Predicts Phase Transition in 

[Ca^li 

Based on prior approximations of heterogeneity in cellular 
excitability and coupling. Boolean networks of connectivity were 
simulated to predict how general multicellular electrical activity 
depends on the relative excitability of the constituent cell 
population and the coupling between individual cells [14,20]. 
Nodes within a cubic lattice had a probability P^xc of being active, 
and adjacent nodes were functionally coupled with a 'coupling 
probability' p (figure lA), and resultant clusters of coupled nodes 
were identified. 'Inexcitable' P-cells can suppress activity in 
excitable P-ceUs via electrical coupling [17], with <30% inexcit- 
able cells necessary for this suppression. To simulate this, a logic 
rule was used for each cluster of coupled nodes within a given 
lattice, where greater than a threshold percentage of inexcitable 
cells (Sp) can suppress activity in all other cells in its coupled 
network. Simulations of the resultant average network activity 
were run with varying values of Pe-xc, Sp, and p to represent 
dilfering cellular excitabUities and electrical coupling (figure IB, SI 
figure SI). 

An increase in electrical activity is predicted 3.S Pfixc 

increased; 

however the functional form is highly dependent onp (figure IC). 
In the absence of coupling (p = 0), a trivial linear response is 
obtained where P^^c represents the level of electrical activity. With 
increasing p, the activity becomes increasingly non-linear as a 
function of P^xc- For higher values oip (0.3 to 1) a sharp transition 
between active and quiescent behavior is observed, representing a 
phase transition with emerging critical behavior. These higher 
values of^ lead to network-spanning coupling (figure SI, S2), and 
as such the 'rule' governing suppression acts over the whole 
network. For low values of ^ (0 to 0.2), the network is composed of 
coupled 'clusters' (figure S2), and the simulation is close to linear 
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Figure 1. Boolean network model predictions. A) Schematic 
representation of the network model with limited connectivity. Note 
larger connected clusters have a higher probability of containing 
inexcitable cells. B) Example false-color maps displaying probability of 
activity, generated from a simulated network with p = 0.30 at Pe^c = 95% 
(top) and Pexc = 60% (bottom). Note substantially increased likelihood of 
activity with the higher Pexc- Further description for how this was 
generated can be found in figure SI. C) Boolean network model 
predictions for the mean percent active cells as a function of proportion 
of excitable cells (Pexc) for varying coupling probabilities p, with a 
threshold fraction of inexcitable cells Sp = 0.1 5. D) Boolean network 
model predictions for the mean percent active cells as a function of 
coupling probability p, for varying proportion of excitable cells Pexc- 
where Sp = 0.1 5. 

doi:1 0.1 371 /journal.pcbi.1 00381 9.g001 



without a strong transition. This level of p corresponds to 
insufEcient coupKng to span the network, which is similar to the 
critical coupling probability (~0.25) in percolation theory [31]. As 
such for j!)>0.25 there are 3 specific regions of emergent network 
behavior: a small (~10%) decrease for Ppx(:>0.85 ('pre-critical 
state'); a rapid ~75% drop at P^xr^O.SS ('critical state'), then a 
small linear decrease for P,,x(<0.85 ('post-critical state'). The 
critical state emerges when P^xc approaches l-Sp; with the 



sharpness of the transition as well as behavior in the pre- and 
post-critical states being dependent on p. 

Overall, this transition can be understood by considering the 
weU-defmed threshold for activity (Sp) and the network-spanning 
connectivity that occurs above the critical coupling probability (p> 
0.25). For values of ¥i.„^~<{l-Sp) there is a gradual decrease in 
network activity with increasing p (figure ID), representing the 
suppressive effect of couphng. For P,.x(:>(l-S/)), the network 
activity remains high, although a slight drop occurs for low levels 
of coupling. Therefore in a general Boolean network, electrical 
coupling is predicted to lead to critical behavior, where a phase 
transition in the network activity occurs as a function of 
constituent cellular activity. 

Phase Transition in [Ca^^i in Kir6.2'^^^° '^^^^'^' Expressing 
Islets 

To test the Boolean network model predictions, we measured 
intracellular free-calcium activity ([Ca^^];) in islets which had 
defined levels of cellular excitability. Islets were isolated from mice 
with inducible, P-cell specific expression of a mutant ATP- 
insensitive Katp channel subunit (Kir6.2'^'^^'"'''^"'''^) under 
Cre'^^-recombinase control [2 7] . Expression of these over-active 
K^XP channels render P-cells functionally inexcitable, causing an 
absence of insulin release, marked hyperglycemia and diabetes 
[27]. Tamoxifen induction of Cre^'^ controls Kir6.2'-'^^'^'''*^"''^ 
expression levels which can be monitored via GFP co-expression, 
leading to both controllable and quantifiable cellular excitabUities 
(SI figure S3). 

At 20 mM glucose islets show [Ca^"""]; which decreased with 
increasing expression of GFP and therefore Kir6.2'^'^^'^"'*^"'''*^, 
similar to model predictions. This showed critical behavior with 3 
specific regions (figure 2A): For low GFP expression <15% (few 
Kir6.2'-'^^''"'^"''''^ expressing cells), [Ca^^; was active over the 
entire islet with similar behavior to wild-type islets lacking GFP 
(GFP = 0) and Kir6.2t^'^^''''''"'"^^ (figure 2A, BI-II, 'pre-critical' 
behavior). Oscillations were almost fully synchronous in each case 
(not shown). For GFP expression at 10-20% there was a sharp 
drop-off in islet [Ca^^;, where small changes in GFP resulted in 
highly disproportionate changes in [Ca^^]; (fi gure 2A, Bill, 
'critical' behavior). Activity was focused to clustered areas of 
synchronization (not shown). For high GFP expression >25% 
(high number of Kir6.2'^'^^'^"'^"'''^'^ expressing cells), islets showed 
sporadic [Ca"'*']; restricted to increasingly smaller clusters (Figur- 
e 2DIV, 'post-critical behavior'). Although islets with <20% GFP 
has similar overall activity compared to wild-type islets (0% GFP), 
there was a marked reduction in the plateau fraction in the <20% 
GFP group (15^0%) compared to the wild-type group (50-70%); 
indicating that even small numbers of inactive cells impacts global 
behavior. In comparison, islets with high levels of GFP (>25%) 
had a low plateau fraction (15-20%) in those cells that were active. 
GFP'*' cells showed similar activity to GFP~ cells albeit with a small 
increase in activity, likely due to a few inactive non-fi-ceUs 
included in the GFP~ analysis (SI figure S4). 

Comparison of experimental data to the Boolean model can be 
seen for a number of values ot p (figure 2C) and Sp (figure 2D). 
Varying p (gap junction coupling) matches the sharpness of the 
transition, whereas varying Sp (number of inexcitable cells 
required to suppress activity) matches the position of the transition. 
A^ = 0.3 (95% CI: 0.280-0.311) and Sjf) = 0.135-0.15 best fits the 
experimental [Ca^""^; data (figure 2A,E). The distribution for fitted 
p was relatively broad but for Sp was well defined (figure 2E). 
These values of p are similar to those found in prior studies 
examining the synchronization of [Ca^"*^; oscillations, which 
indicated a limited level of functional coupling in the islet 



PLOS Computational Biology | www.ploscompbiol.org 



3 



September 2014 | Volume 10 | Issue 9 | el 00381 9 



Phase Transitions in Islet Excitability 



IV 



B 



100 
80 
60 
40 

20 



^0.3 
(D 
= 0.2 

a; 0.1 





100 


JO 


80 


aj 




O 


60 


0) 




.> 






40 


< 








20 




0 



■ p = 0.30 

■ Sp = 0.15 
-95% CI 




0.4 0.6 
P (1-GFP) 




100 



"0 
O 

(U 



80 



60 



t3 40 
< 

20 




-Sp = 


0.05 


-Sp = 


0.10 


-Sp = 


0.15 


-Sp = 


0.20 


-Sp = 


0.25 




0.2 0.4 0.6 O.i 

P (1-GFP) 



0.2 0.4 0.6 O.i 
P (1-GFP) 



IVI 





0.26 0.28 0.3 
P 



60 120 180 240 300 
Time (Seconds) 



Figure 2. Experimental data sKiowing how Boolean network model describes phase transitions in islet [Ca^"^]!. A) Percent cells showing 
[Ca^*]i elevations as a function of number of excitable cells, as determined by lack of GFP and thus Kir6.2'*'^^°''^^^^°' expression (i.e. Pexc = 1-%GFP), 
together with Boolean network model fit. Filled squares indicate experimental data, solid line represents mean of simulations that best fit data with 
p = 0.30 and Sp = 0.15 (x^=1.38), dashed lines represents 95% confidence intervals of the simulation fit. B) Representative [Ca^*]i data for islets 
indicated in A, from regions of wild-type (I), 'pre-critical' (II), 'critical' (III) and 'post-critical' (IV) levels of Pe^c- Left: Areas of activity are highlighted in red 
and scale bars represent 50 |im. Right: Representative time-courses of normalized FuraRed calcium dye fluorescence for cells within each islet, where 
vertical scale bar indicates 20% change in fluorescence. Red time-courses are determined to be active, black time-courses are determined to be 
inactive. See SI for Movies SI, S2, S3, S4 of these data. C) Experimental data with Boolean network simulations for varying connectivity p. D) As in C for 
varying threshold of inactive cells Sp. E) Probability distribution of fitted p (linear scale) and Sp (log scale) parameters to data in A, along with heat 
map of 2D distribution (log scale). 
doi:1 0.1 371 /journal.pcbi.1 00381 9.g002 



()f) = 0.31-0.36) [20]. These values oi Sp are also consistent with 
experimental studies that suggest between 1 and 30% of inactive 
cells can suppress activity in other cells [17,20]. 

Therefore, introducing inexcitable cells into the islet experi- 
mentally generates critical behavior which quantitatively agrees 
with a Boolean network model and predicts the importance of 
electrical coupling in regulating multicellular excitability. 

Phase Transition in Physiological Parameters 

P-cell [Ca""^i drives insulin release to regulate glucose 
homeostasis. Given that the behavior in [Ca^"*]; following varied 
expression of over-active Katp channels, we next tested whether 
this also occurred in downstream physiological parameters. 
Averaged over each mouse, similar [Ca^"*^; was observed in 
wild-type islets lacking GFP and islets with low-level GFP (<20%, 
'pre-critical'), while both were significantly greater than [Ca^"*]; in 
islets with high GFP expression (>20%, 'post-critical') (figure 3A). 
Plasma insulin also showed a similar transition, with pre-critical 
(GFP<20%) plasma insulin being significantly greater than post- 
critical (GFP>20%) plasma insulin (figure 3B). However mice 



lacking GFP did show significantly greater insulin than mice with 
low-level GFP, correlating with the reduced plateau fraction 
observed. Insuhn reduces glucose levels, and as expected pre- 
critical mice (GFP<20%) had normal glucose levels (figure 3C), 
while post-critical mice (GFP>20%) demonstrated elevated 
glucose levels. Glucose-stimulated insulin secretion from isolated 
islets showed similar behavior to that of plasma insulin 
(figures 3D-F), where again islets lacking GFP showed significant- 
ly greater GSIS than islets with low level GFP. 

Therefore insulin dynamics and blood glucose levels follow 
similar behavior as the driving [Ca^"""]; following varied expression 
of Kir6.2'^'^^'^'''''"'^'^, demonstrating a physiological link in the 
critical behavior in [Ca^^]; activity as a function of P^xc,- 

Coupling Dependent Suppression in Kir6.2^'^'^'^' 
Expressing Islets 

The Boolean model accurately predicts the impact of variable 
cellular excitabilities (PexJ on [Ca^'^i suppression at elevated 
glucose through expression of over-active Katp channels 
(Kir6.2'^'^^'^"'^"''^''^). However, the Boolean model also predicts 
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Figure 3. Link between phase transitions in [Ca ]i and 
physiological parameters. A) Percent cells showing [Ca^*]j elevations 
averaged over islets from each Kir6.2"^'^^°''^^^^°'-expressing mouse as a 
function of Pexc (100%-%GFP). Right: IVIean(±s.e.m.) for data binned to 
wild-type, pre- and post-critical ranges as determined by %GFP. 
*indicates significant difference (p<0.0001) between data as indicated. 
B) Plasma insulin levels from each mouse as a function of Pexc- Right: 
IVIean(±s.e.m.) for data binned as in A. 'indicates significant difference 
(p<0.05) between data as indicated. C) Time-averaged blood glucose 
levels from each mouse as a function of P^xc- Right: Mean(±s.e.m.) for 
data binned as in A. 'indicates significant difference (p<0.0001) between 
data as indicated. D) Insulin secretion from isolated islets at 20 mM 
glucose (left) and 2 mM glucose (right), averaged over each mouse as a 
function of P^xc- E) Mean(±s.e.m.) for data in D binned as in A. 'indicates 
significant difference (p<0.01) between data, 'ns' indicates no significant 
difference (p>0.05) as indicated. F) Mean(±s.e.m.) of islet insulin content 
averaged over each mouse. Solid lines in A,B,D represent Boolean model 
fit (p = 0.3) from experimental data in Figure 2C. 
doi:1 0.1 371 /journal.pcbi.1 00381 9.g003 



how [Ca"^i .suppression varies as a function of gap junction 
coupling (p) (Figure ID). To test tliis, we measured [Ca^"*]; in islets 
from mice with P-cell specific mosaic expression of an inactive 
Katp channel subunit (Kir6.2[^^ This was combined with a 
knockout of Cx36, yielding 100% (Cx36'^^\ 50% (CxSe""^") or 
0% (Cx36 '' ) gap junction coupling, as well application of the 
gap junction inhibitor 1 8-a-glycyrrhetinic acid [10]. Expression of 
inactive Katp channels render P-cells constitutively (glucose- 
independent) active, yet islets which have a majority (but not all) of 
their cells expressing inactive Katp channels show glucose- 
dependent electrical activity similar to wild-type islets [17]. GFP 
co-expression indicates ~70% of P-ceUs express inactive Katp 
channels (SI figure S3) such that Pexc = 0.7. With increasing gap 
junction coupling [Ca""^; progressively decreased until residual 
activity was observed at full coupling, similar to that in the post- 
critical state upon Kir6.2'-'^^^'''^"'^^ expression. There was strong 
agreement between experimental measurements and the Boolean 
Network model, with ap at normal (0x36''"''''", 100%) gap junction 
coupling of 0.38 (95% CI: 0.372-0.394) and a suppression 
threshold Sp = 0.15 giving the best fit (figure 4A). This is similar 
top, Sp derived in the first experimental system (figure 2). Varying 
Sp affects the gap junction dependence in [Ca^^;, with little effect 
between 0.05-0.2, but strong divergence above 0.2 (figure 4B). 
Therefore the Boolean network model can accurately predict 
behavior in a different experimental model with defined levels of 
cellular excitability (Pexi) and gap junction coupling Ip). 

Phase Transition in a Dynamical Model of p-Cell 
Networks 

The Boolean network model accurately describes how [Ca^^; 
critically depends on cellular excitability and coupling. Neverthe- 
less it is a static framework of a dynamical system and does not 
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Figure 4. Boolean network model describes [Ca^"^]! suppression 
as a function of coupling conductance. A) Percent cells showing 
[Ca^*]i elevations in islets from Kir6.2"^'°''^'-expressing mice as a function 
of gap junction conductance, together with Boolean network model fit. 
Filled squares indicate mean(±s.e.m.) experimental data, solid line 
represents mean of simulations that best fit data for wild-type coupling 
value po = 0.38 and Sp = 0.1 5 (x^ = 0.416), dashed lines represents 95% 
confidence intervals of simulation fit. Gap junction conductance for 
each data point was normalized to the wild-type conductance and 
scaled by the fitted pg. B) Mean(±s.e.m.) experimental data with 
Boolean network simulations for varying threshold of inactive cells Sp. 
doi:1 0.1 371/Journal.pcbi.1 00381 9.g004 
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take into account limit-cycle behavior. To investigate whether 
similar behavior exists in a coupled dynamical oscillator model of 
the islet, we generated a multi-cellular version of a recent P-ceU 
model which includes a comprehensive description of P-ceU 
electrophysiology [30]. Our model also included a more realistic 
quasi-spherical architecture, heterogeneity in gap junction cou- 
phng [14,32], and heterogeneity in endogenous cellular activity 
[11,14]. To model K^-j-p-overactivity resulting from 
Kir6.2'-'^^"''*^"'''^ expression, a defined fraction of cells with 
reduced ATP-inhibition of Katp activity was introduced to render 
them inexcitable. 

As with the Boolean network model and experimentally 
measured [Ca^'"'];, a clear phase transition was observed at 
20 mM glucose in the coupled oscillator model with ~15% 
KATP-overactivity (figure 5A). Again critical behavior manifested 
in three regimes. Simulated islets without Katp over-activity 
showed [Ca""^; dynamics closely matching previously published 
models (figure 5BI) [30]. Simulated islets with low KATP-overac- 
tivity (<15%) showed a linear decrease in activity with a reduced 
plateau fraction as experimentally observed (figure 5A, BII, 'pre- 
critical' behavior), while maintaining near-full synchronization. 
Simulated islets with KATP-overactivity at 10-30% again showed a 
sharp drop-off in [Ca^"*^;, with small changes in KAxp-overactivity 
leading to highly disproportionate changes in [Ca^-*]; (figure 5A, 
Bill). Simulated islets with high KATp-overactivity (>30%) 
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Figure 5. Coupled dynamical oscillator model describes 
experimental islet phase transitions. A) Percent cells showing 
[Ca^*]i elevations in simulated islets as a function effraction of excitable 
cells (Pexc). as set by the % cells lacking ATP-insensitivity. Solid line 
represents mean of simulation results generated from 5 random 
number seeds, dashed lines represents 95% confidence intervals of 
simulations. B) Representative simulated [Ca^*], time-courses for 
parameters indicated in A, from regions of wild-type (I), 'pre-critical' 
(II), 'critical' (III) and 'post-critical' (IV) behavior, as in figure 2. Vertical 
scale bar indicates 20% change in simulated [Ca'^*]j. Red time-courses 
are determined to be active, black time-courses are determined to be 
inactive. See SI for Movies S5, S6, S7, S8 of these data. C) Percent cells 
showing [Ca^*]j elevations in simulated islet as a function of number of 
excitable cells (Pgxc) for varying mean gap junction conductance values. 
Filled squares indicate experimental data from Kir6.2'^'^^°''*'*'^°'- 
expressing islets in figure 2. 
doi:1 0.1 371/journal.pcbi.1 00381 9.g005 



showed only sporadic low level [Ca"^; (Figure 5 A, BIV, 'post- 
critical' behavior). A physiological mean gap junction conductance 
of 120 pS [14,32] was found to best describe experimental data 
(figure 5A,C). The sharpness and position of the phase transition 
was highly dependent on the mean coupling conductance, with 
increasing conductance leading to a sharper transition occurring at 
lower KATP-overactivity (figure 5C). 

The islet is commonly modeled as a cubic lattice or other 
regular geometry [14,21,33]. A spherical islet-like structure which 
has a heterogeneous number of cell-cell connections 
(meaii,SD = 5.3,1.7) generated a less-sharp transition compared 
to a regular cubic geometry, and this better matched experimental 
data (SI figure S5). Similarly, a heterogeneous level of coupling 
conductance generated a less-sharp transition (SI figure S6A). This 
indicates the importance of coupling heterogeneity, in terms of 
connection geometry, connection number and connection 
strength. The endogenous heterogeneity in cellular activity did 
not significantly impact the phase-transition indicating the 
dominating effect of Krr6.2'^'^^'^'''^"*''*^ expression (SI figure 
S6B). A similar phase-transition was also observed for simulations 
run at 1 1 mM glucose (not shown). 

Therefore critical behavior also emerges in a dynamic coupled 
fi-cell oscillator model with quantitative agreement with experi- 
mental measurements and a static Boolean network model. 

Phase Transition Resulting from Endogenous 
Heterogeneity and Coupling 

We have examined how the coupling between heterogeneous 
cells leads to critical behavior by introducing defined mutant 
populations of inexcitable ceUs (Kir6.2[^^'''''''"'-^'2]-) ^j. excitable 
cells (Kir6.2'^'^^^). However, endogenous P-cells are themselves 
highly heterogeneous under physiological ranges of glucose, 
showing a continuum of excitabUities rather than being constitu- 
tively excitable /inexcitable [9-11]. 

To examine how gap junction coupling leads to critical 
behavior in the presence of endogenous heterogeneity, we applied 
a 'ramp' of increasing diazoxide concentrations to uniformly 
promote Katp channel opening. At 1 1 mM glucose, [Ca^'']; in 
wild-type islets at 0 |J,M and 50 |J.M diazoxide was similar, but at 
100 |J.M there was a rapid ~60% drop (figure 6A), where only a 
few remaining cells were active (figure 6B). Similar low-level 
[Ca^+j i was observed at 250 |iM. In islets from mice lacking Cx36 
gap junction coupling, similar [Ca^"^; was observed to wild-type 
islets at 0 |XM diazoxide, albeit with no synchronization. Upon 
increasing diazoxide, a more gradual decrease in [Ca^"*]; was 
observed, with less [Ca^"^; observed at 50 |J.M diaozrxde but more 
[Ca^'^i remained at 100 |lM diazoxide (figure 6A,B). 

These data were also well described using the coupled dynamic 
oscillator model. In the presence of endogenous heterogeneity at 
1 1 mM glucose, a uniform reduction in ATP-sensitive Katp 
inhibition led to a clear phase transition in islet [Ca^""^; in the 
presence of normal coupling (120 pS) (figure 6C,D). However, in 
the absence of coupling a more gradual change occurred in good 
agreement with experimental measurements (figure 6C,D); where 
[Ca^-^i was elevated in the absence of coupling over a certain 
range of uniform Katp inhibition. As such, <50 [iM diazoxide lies 
in the 'pre-critical' regime, > 1 00 |J.M diazoxide lies in the 'post- 
critical' regime, and the transition lies at 50—100 |J,M. 

In experiments with mutant Katp subunit expression, cells were 
considered 'inexcitable' if they showed GFP and 
Krr6.2'-'^^''"'*^"'''*^ expression. In this case of endogenous hetero- 
geneity, for a given concentration of diazoxide, we can consider a 
cell is 'inexcitable' if it is quiescent in the absence of electrical 
coupling. By plotting activity in the presence of coupling 
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Figure 6. Coupling-dependent phase transitions upon uniform Katp inhibition resulting from endogenous heterogeneity. A) 

Mean(±s.e.m.) percent cells showing [Ca^*]| elevations at 2 mlVI glucose and at 11 mM glucose with 0, 50, 100, 250 |iM of the K^jp activator 



diazoxide, for wild type islets (Cx36 , solid squares) and islets lacking gap junction coupling (Cx36 , empty circles). 



indicate significant 



difference (p<0.05, p<0.005) between activity in Cx36 and Cx36" islets at treatments indicated. B) Representative [Ca time-courses for 
Cx36*'* and Cx36^'^ islets at 1 1 mM glucose with 0 |iM or 100 |iM diazoxide. Red time-courses are determined to be active, black time-courses are 
determined to be inactive. Vertical scale bar indicates 20% change in fluorescence. C) Percent cells showing [Ca^*], elevations in simulated islet as a 
function of a uniform increase in the fraction of ATP-insensitivity of Katp channel activation (a) across cells of the islet. Mean simulation data is 
presented for zero gap junction conductance (0 pS) and wild-type gap junction conductance (120 pS). D) Representative simulated [Ca^*], time- 
courses for wild-type gap junction conductance (120 pS) and zero gap junction conductance upon indicated levels of uniform Katp activation. Red 
time-courses are determined to be active, black time-courses are determined to be inactive. Vertical scale bar indicates 20% change in simulated 
[Ca^^, 

doi:1 0.1 371 /journal.pcbi.1 00381 9.g006 



(representing the resultant activity) against activity in the absence 
of coupling (representing intrinsic cellular excitability) similar 
phase-transitions are apparent; with quantitative agreement 
between experimental data, dynamic coupled oscillator model 
and static network model (figure 7A-C). The phase transition in 
the dynamic coupled oscillator model was dependent on how 



heterogeneity was generated, where heterogeneity in multiple 
factors rather than any one factor was required for agreement with 
experimental data (SI figure S7). Therefore critical behavior can 
occur more generally from the coupling between heterogeneous 
cellular populations within the islet, as exemplified here experi- 
mentally and theoretically. 



ABC 




Figure 7. Phase transitions in endogenous p-cell network activity, as shown by the activity in a fully-coupled islet system as a 
function of the activity in the uncoupled islet system; where the latter represents the intrinsic excitability of the constituent cells. A) 

Experimentally measured transition from global activity to quiescence in wild-type islets treated with varying diazoxide concentrations, showing 
phase transition in activity as constituent cellular activity is reduced B) Simulated transition from global activity to quiescence upon normal gap 
Junction conductance as Katp is uniformly activated across the islet in the dynamical oscillator model. C) Modelled transition from activity to 
quiescence within the Boolean lattice resistor network model as Pexc is reduced, for p = 0.3 and Sp = 0.5. Note in all cases; for islets lacking gap 
Junction coupling, with zero gap junction conductance and for p = 0, the transition is trivially linear (blue dashed). 
doi:1 0.1 371/Journal.pcbi.1 00381 9.g007 
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Discussion 

The islet of Langerhans shows unique functional properties that 
result from the underlying network interactions between constit- 
uent cells. One important property is that |3-cells within the islet 
show globally quiescent behavior at lower levels of glucose despite 
showing a heterogeneous range of glucose sensitivities when in 
isolation. A proposed mechanism underlying this behavior is that 
at a given glucose stimulation, inactive cells suppress cells that 
otherwise would be active, via gap junction coupling. We applied 
predictive mathematical models to quantify this behavior and 
determined the relative role of Katp channel activity (controlling 
cellular excitability) and gap junction activity (controlling cellular 
coupling) in shaping this islet response. We then experimentally 
verify predicted behavior using two independent experimental 
models. 

Critical Behavior in Islets and Theoretical Models 

In hne with previous work describing coupled electrical 
dynamics, we showed that the structure and function of the islet 
cellular network can be described through principles of network 
theorV' [20,21]. Both the Boolean network and dynamic oscillator 
m<)d('ls predict the emergent behavior upon coupling between a 
heterogeneous cell population. The islet rapidly transitions 
between globally coordinated active and inactive states upon 
disproportionally small changes in the excitability of the constit- 
uent cells as they approach a critical 'threshold' excitability. This 
occurs under both conditions of P-cell heterogeneity we examined: 
the imposed bimodal P-cell populations achieved through 
expression of Kir6.2t'^^°'^"'^^ or Kir6.2^^^ mutations; and 
endogenous P-ceU heterogeneity with diazoxide activation of 
Katp- 

The Boolean model reveals that there is an imbalance in the 
ability of excitable and inexcitable cells to respectively propagate 
stimulation or suppression. A low Sp in the model indicates a 
preference for excitable cells to be suppressed by inexcitable cells. 
This describes how gain-of-function Kir6.2'^'^^'''^'^^^ expressing 
cells (which are glucose-unresponsive) suppress activity in coupled 
normal cells at high glucose, and how loss-of-function Kir6.2'^'^^^ 
expressing cells are suppressed by normal cells at low glucose 
(figure 2, 4). The role of p (gap junction coupling) determines the 
spatial extent over which suppression occurs. As shown in figure 
SI and S2, a low^ results in coupled behavior restricted to a few 
cells and therefore inactive cells are unlikely to couple to many 
active cells and mediate suppression. When p exceeds the critical 
coupling probability (~0.25) then coupling spans the whole 
network and inactive cells can couple to and suppress most active 
cells in the network. The sharp transition that emerges upon p> 
0.25 can be understood by considering that the threshold for 
activity (Sp) is well defined with a sharp cutoff for the V^^^ which 
determines whether the cluster is active or inactive. The 
agreement with experimental data indicates that there is little 
variability between cells in this threshold for suppression, as also 
supported by the distributions of fitted Sp (figure 2E). 

While the coupled dynamic oscillator model also predicts and 
describes the phase transitions present, the Boolean model 
describes the essential features that govern multicellular regulation 
of islet excitability. Results suggest that the islet may fundamen- 
tally behave in a binary fashion in terms of gap junction coupling 
and KATP-regulated excitability. Given the proportion of cells that 
intrinsically (i.e. in the absence of coupling) show activity at a given 
glucose stimulation and the level of coupling, the overall response 
of the islet can be approximated through this reductionist model. 
Of course dynamical features are missing from the Boolean model 



which is only described by the coupled dynamic oscillator model: 
including the altered oscillator)- characteristics in the pre-critical 
state. The low p (0.30-0.38) required for the Boolean network 
model to quantitatively describe experimental data points to 
incomplete coupling present; and this can explain the residual 
activity in the post-critical state (figure IC). Recent studies of 
coordinated [Ca^^; oscillations and waves in the islet have 
indicated a 'backbone' of a few strong connections dominate 
coupled [Ca^"']i dynamics, which is equivalent to a similarly low p 
[19,20]. The ability of the coupled dynamic oscillator model to 
also describe the transition between globally active and inactive 
states, suggests that the dynamics of the islet may behave 
according to general principles of coupled dynamical systems. 
Further work is needed to examine this critical behavior in more 
detail, including power law scaling and its dependence on network 
parameters and cellular properties. 

The phase-transition behavior can also be explained through a 
mean-field theory analogy (SI figure S8). Cells expressing the 
mutant Kir6.2'^'^^^°'^"'^'^ are intrinsically inexcitable (figure 
S8A). In the 'pre-critical' regime the number of these cells is 
below a critical threshold and insufficient to suppress glucose- 
stimulated activity via coupling; therefore all cells are recruited to 
ele\'ate [Ca^^;. When the number of these inactive cells 
approaches the critical threshold (Sp = 0.15 for Kir6.2'-'^^'"''*^"'^^ 
expression) critical behavior emerges and coupling mediates 
suppression of other active cells. In normal islets treated with 
diazoxide, endogenous P-ceU heterogeneity leads to variable 
intrinsic excitabilities and we expect diazoxide renders cells less 
glucose sensitive to be inexcitable (figure S8B). In the absence of 
couphng these are observed to be inactive (figure 6). Low 
concentrations of diazoxide (<50 JtM) render only a few cells 
inexcitable, which is below the critical threshold (Sp~0.5) and 
insufficient to suppress [Ca^""];. At higher concentrations of 
diazoxide (> 1 00 |J,M) more cells are rendered inexcitable, and 
when this number exceeds the critical threshold, coupling 
mediates suppression of other normally excitable cells. We predict 
that observed glucose-dependent activity and the coupling 
dependence can also be explained in this way (see below). 

The Sp for endogenous heterogeneity is higher than that for an 
imposed biomodal distributions (i.e. diazoxide treatment versus 
Kir6.2'^'^^^'''^'^^^ or Kir6.2'^'^^^ expression) suggesting a more 
even balance between the ability of excitable and inexcitable cells 
to respectively propagate stimulation or suppression in wild-type 
islets (figure 6, 7). This balance may arise from the difierent 
distribution of heterogeneity present, but a phase-transition stiU 
emerges in the presence of coupling indicating a more general 
regulation of multicellular excitability. 

Therefore through limited coupling of heterogeneous popula- 
tions of cells, critical behavior emerges in the islet dynamical 
system where large changes in activity result from small changes in 
the constituent cellular excitabilities. 

Applicability to Physiology and Diabetes 

Gap junctions impact islet behavior in two main ways. At high 
glucose (Katp channel-closure), gap junctions coordinate oscilla- 
tory dynamics of membrane depolarization and [Ca^^i to 
generate a robust pulsatile insulin secretion [13,14,32]. A number 
of recent studies have examined this aspect, including multicellular 
modelling and quantitative analyses [14,18,19,21]. Equally 
important however is that at lower glucose (Katp channel- 
opening), gap junctions mediate a suppression of membrane 
depolarization, [Ca^""^;, and insulin secretion [10,13,17]. The 
mechanisms involved in mediating suppression are not well 
characterized, and several experimental perturbations have 
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yielded unexpected results or have not been well described 
theoretically [17,34,35]. Here, we were able to quantitatively 
describe suppressive behavior resulting from coupling, which 
yields a more complete understanding for how the islet functions 
under conditions of K^i p channel opening. 

At 6—7 mM glucose, the islet sharply transition between global 
quiescence and globally synchronized [Ca^"^i oscillations. In the 
absence of coupling, the progressive elevation in the number of 
cells showing [Ca^"^]; elevations is gradual [10]. This follows the 
same behavior as variable Kir6.2'^'^'^^''''^'^^'^ expression and 
diazoxide concentration (figures 2,6). At <6 mM glucose, global 
suppression is equivalent to >15% Kir6.2'^^^'''''''"'"^''^ expression 
or > 100 (iM diazoxide; whereas at >7 mM glucose global activity 
is equivalent to <15% Kir6.2^^'^^°''''^*'2] expression or <50 \lM 
diazoxide. The 6-7 mM glucose transition is therefore equivalent 
to behavior at ~15% Kir6.2'-'^^'^"''^'^''''^ expression or ,50- 
100 jxM diazoxide. As such, we propose results from the Boolean 
model, as illustrated by the mean-field theory analogy, have 
greater implications by describing physiological glucose-dependent 
islet electrical activity (figure S8C). Coupling heterogeneity and 
islet architecture lead to variability in the number and strength of 
connections, impacting the phase transition. These factors may 
therefore play a role in shaping the physiological regulation of 
glucose-stimulated [Ca^^; and insulin secretion (Figures S5,S6). 

At 1 1 mM glucose, heterogeneity leads to a small population of 
P-cells (<10%) remaining inactive in the absence of coupling [10]. 
In the presence of couphng there is global activity with a lower 
plateau fraction compared to higher glucose levels (e.g. 20 mM) 
[36]. This matches the behavior at 5-10% Kir6.2[^^'"'''"''-^'2] 
expression or 50 |J,M diazoxide in the respective absence and 
presence of couphng. Therefore an alternative view for how 
oscillatory dynamics are shaped at an islet-wide level is that less- 
active cells within the fi-cell network have a modulatory effect on 
overall oscillation waveform, rather than oscillaticms being shaped 
by purc'ly intrinsic properties of the fi-cells. Importantly, the 
reduced plateau fraction of [Ca "*]; bursts at 5-10% 
Kir6.2^'^^ ,Ki85Q] expression correlates with a significant decrease 
in insulin secretion (figure 3). A decrease in burst duration has 
previously been suggested to reduce insulin release [37], as 
supported by these results. Thus subtle alterations in the balance of 
constituent cell excitabHities have a strong physiological effect on 
islet function. 

Our results also have implications for neonatal diabetes meUitus 
(NDM), where the majority of cases result from mutations to 
Kir6.2 or SURl Kaip channel subunits [8,38]. Kir6.2f^^''"''^"'^'2^ 
expression models this disease [27]. Our results show that NDM 
mutations gives rise to a disproportionate suppression in [Ca^"']i 
and insulin release, thereby causing diabetes due to the critical 
behavior that emerges from coupling and network dynamics. This 
also explains how an absence of coupling elevates [Ca^^j and 
insulin release (figures 1, 6) to prevent the progression of diabetes. 
This was experimentally demonstrated in a recent study [26], and 
the rescue of diabetes can only be understood mechanistically at 
the multicellular level. Clearly in human diabetes, mutations are 
not expressed mosaically. However, the diazoxide results which 
depend on a continuum of heterogeneity (figure 6) demonstrate 
that critical behavior exacerbates NDM upon uniform K\ip 
channels overactivity. Other monogenic diabetes causing muta- 
tions that affect P-ceU excitability, such as Glucokinase [6], may 
also have similar effects on islet excitability and lend themselves to 
analysis by the Boolean model and coupled oscillator model. 

Mutations causing NDM are functionally equivalent to >15% 
Kir6.2'^'^^''"'^"'^^ or > 100 |iM diazoxide, effectively residing in a 
post-critical state suppressing global [Ca^"*];. There exists a 



spectrum of Katp channel mutations linked to diabetes, where 
weaker mutations to Kir6.2 and SURl elevate the risk of type2 
diabetes [39-41]. These mutations likely have a more subtle effect 
on islet excitabilit}' and as a consequence we predict that islets 
residing in the pre-critical state (<15% Kir6.2L^*^'*"'^"'''21 or < 
50 |J,M diazoxide) would still be susceptible to diabetes following 
metabolic stress. Further, while gap junction reduction recovers 
insuUn release and glucose control in the post-critical regime (i.e. 
NDM), we predict a gap junction increase would be beneficial in 
the pre-threshold regime (i.e. type2 diabetes). 

Converse to this, results from Kir6.2'^'^^-' islets show how 
critical behavior provides the islet with a resilience to over- 
excitable P-cells. Given the ~85% threshold of excitable cells 
required to elevate [Ca^"*]!, with ~70% over-excitable Kir6.2^" 
^-expressing cells, many of the ~30% normal P-ceUs would 
also need to be active (e.g. >50% at ~5.5 mM glucose [10]). This 
explains only the minor shift in glucose-stimulated [Ca^^; that 
occurs following Katp inactivity and highlights the role electrical 
coupling plays in protecting islets against hyper-excitabUity [17]. 

Therefore we describe the emergence of critical behavior linking 
multiple levels including molecular and ( (-Uular behavior, multi- 
cellular behavior, in-vivo physiology, disease and treatment. 

Implications for General Dynamical Systems 

This study also has implications for general understanding of 
physiological systems composed of coupled dynamic units. 
Previous theoretical studies have shown how the introduction of 
non-oscillatory elements above a critical level in a generalized 
coupled oscillator system can lead to cessation of global oscillations 
with phase transitions [42,43]. Our study experimentally and 
theoretically demonstrates this in a disease relevant system. 
Further, prior studies theoretically demonstrated that the fraction 
of excitable elements (i.e. P,.x< ) and coupling strength (i.e. p) exist in 
a phase plane where increased coupling decreases the number of 
inactive elements required for suppression [43,44]. We demon- 
strated this experimentally and theoretically, with 
Kir6.2'^'^'^^'''^'^^^ and diazoxide-induced suppression. While 
strong coupling promotes robust synchronization, it will increase 
suppression from non-oscillatory inexcitable units. Given a small 
population of inactive P-cells exists as a result of cellular 
heterogeneity, these generalized theoretical studies imply inappro- 
priate elevations in coupling would be deleterious, by reducing 
glucose-stimulated [Ca^"']i. As such the level of coupling is likely at 
an optimal level to balance global synchronization and suppres- 
sion. 

The strong link between dynamical P-ceU networks and 
generalized coupled oscillators implies similar behavior can be 
expected in other physiological systems. In the heart, electrical 
activity is initiated by pacemaker cells and propagates to excite 
contractile myocytes. In culture, non-excitable fibroblasts propor- 
tionally reduced cardiomyocyte wave propagation bursts frequen- 
cy with Cx43 dependence [45]. No activity was reported for > 
30% fibroblast penetrance and modulation of action potential 
frequency occurred at <30% fibroblast penetrance (implying 
Sp = 0.3). While the mechanisms of coupling-dependent suppres- 
sion are very different compared to our study, global responses are 
similar implying similar governing principles. Similarly, pacemak- 
er cells exhibit dominance over myocytes at an optimal gap 
junction conductance [46]; where high couphng leads to 
arrhythmias and low coupling leads to poor synchronization 
[47]. Neurons also display intrinsic oscillatory behavior, and the 
effects of coupling and presence of inhibitory and excitatory 
neurons on synchronization and phase modulation is an active 
area of research [48-53]. Critical dynamics have been described 
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theoretically to emerge from excitatory and inhibitory units in 
neuronal networks [54], and a computational study which 
introduced 'contrarian elements' into neural coupled oscillator 
networks found that a similar threshold of 15% suppressed global 
dynamics [49]. 

We also anticipate that critical behavior resulting from coupling 
of heterogeneous units may be considered a general regulatory 
mechanism. Many systems respond to a stimulus by transitioning 
between inactive and active states (e.g. contractile, hormone- 
secretory). Our study implies that constituent cellular units need 
not themselves have a uniform or robust response to generate a 
robust multicellular response. Rather, a robust response can 
emerge from coupling a heterogeneous collection of cells; where 
coupling and architecture need only have sufficient strength and 
connection number on average. This makes the overall system 
robust against noise and variability, and loosens the requirement 
for tight regulatory mechanisms within the constituent cells. 
Similarly, given a constant stimulus, a robust transition between 
globally active and inactive states could be achieved by remodeling 
connectivity with a small number of inhibitory units. For example, 
down-regulating connections from ~20% to ~5% inhibitory cells 
would transition the system from inactive to active requiring 
minimal system remodeling. This also suggests how inappropriate 
changes in coupling or constituent cells may lead to global non- 
responsiveness and disease. We speculate these principles may 
apply to other neuroendocrine cell systems such as GH-cells or the 
adrenal medulla, where functional remodeling elevates hormone 
secretion upon physiological stimuli [18,5,5]. Indeed many of these 
principles have been linked with GnRH neuron function during 
development [56,57]. Therefore we suggest that new and robust 
functionalities can be generated at the multicellular level from the 
coupling of non-robust constituent cell function, requiring minimal 
system resources compared to the requirements were cells to act 
autonomously. 

Conclusions 

Any living system cannot avoid deterioration through mutation 
or other pathological insult. This study experimentally and 
theoretically demonstrates that if the fraction of inactive elements 
exceeds a coupling-dependc'ut threshold, the global activity of the 
system can be abolished. In the case of the islet this explains how 
inactive cells can suppress the activity of other cells, thereby 
preventing the secretion of insulin. In the case of Katp mutations, 
this quantifies the threshold of inexcitable cells required for 
pathogenic symptoms and explains how coupling can eliminate 
the emergence of diabetes oi exacerbate it. Overall, this gives a 
new understanding for how emcrg<;nt properties of the islet as a fS- 
ceU network are generated; as well as for understanding islet 
dysfunction in diabetes and novel ways to overcome dysfunction. 
More broadly, this generates insight into emergent behavior of 
multicellular systems in general. 

Methods 

Ethics Statement 

All experiments were performed in compliance with the relevant 

laws and institutional guidelines, and were approved by the 
University of Colorado Institutional Biosafety Committee (IBC) 
and Institutional Animal Care and Use Committee (lACUC). 

Animal Care 

The generation of Rosa26-Kir6.2f^^^"'''"''^'2] ('gain-of-func- 
tion' Katp subunit with GFP co-expression), Pdx-Cre^^ (fS-cell 
specific inducible Cre), Kir6.2'^^^^ ('loss-of-function' Katp 



subunit with GFP tag), and Cx36 (Connexin36 global 
knockout) have been described previously [27,28,58,59]. Expres- 
sion of variable Kir6.2'^^^^''''^' was achieved in P-cells by 

crossing Rosa26-Kir6.2'"'^^''"''^'^''''^ and Pdx-Cre^'*^ mice, and 
inducing Kir6.2'^'^^^"''^"'''''^ expression in 8-16 week old mice by 
1-5 daily doses of tamoxifen (50 mg/g body-weight). Littermates 
lacking Rosa26-Kir6.2'^'^^''°'*'"'^*^ and/or Pdx-Cre^^ were used 
as controls. 

In-Vivo Measurements 

Blood glucose was measured daily and averaged over day 27-29 
post tamoxifen induction using a glucometer (Ascensia Contour, 
Bayer). Plasma insulin was measured at day 29 from blood samples 
centrifuged for 15 minutes at 13,900RCF, and assayed using 
mouse ultrasensitive insulin ELISA (Alpco). 

Islet Isolation and Insulin Secretion 

Islets were isolated by coUagenase injection into the pancreas 
through the pancreatic duct; the pancreas was harvested and 
digested, and islets were handpicked [23]. Islets were maintained 
in RPMI medium (Invitrogen) supplemented with 10% FBS, 
1 1 mM glucose, 100 U/ml penicillin, 100 Hg/ml streptomycin, at 
37°C under humidified 5% CO2 for 24—48 hours prior to study. 
For insulin secretion measurements, islets (5 /column, duplicates) 
were pre-in[:ubated in Krebs-Ringer buffer (128.8 mM NaCl, 
5 mM NaHCOs, 5.8 mM KCl, 1.2 mM KH2PO4, 2..5mM 
CaCla, 1.2 mM MgSO*, 10 mM HEPES, 0.1% BSA, pH 7.4) 
plus 2 mM glucose; then incubated for 60 minutes in Krebs- 
Ringer buffer plus 20 mM glucose. After incubation, the medium 
was sampled and insulin concentration assayed using mouse 
ultrasensitive insulin ELISA. To estimate insulin content, islets 
were lysed in 1% TritonX-100 and frozen at — 20°C overnight. 

Calcium Imaging 

To measure [Ca^^; dynamics, isolated islets were loaded with 
4 |j.M FuraRed-AM (Invitrogen) in imaging medium (125 mM 
NaCl, 5.7 mM KCl, 2.5 mM CaCla, 1.2 mM MgClg, 10 mM 
Hepes, 2 mM glucose, and 0.1% BSA, pH 7.4) at room 
temperature for 90-120 minutes and held in polymdimethylsilox- 
ane (PDMS) microfluidic devices [17]. FuraRed fluorescence was 
imaged on a spinning-disk confocal microscope (Marianas, 31) 
with a 40x 1.3NA Plan-NEOFluar oil-immersion objective (Zeiss) 
maintained at 37°C. Images were acquired at 1 frame/ sec using a 
488 nm diode laser for excitation and a 580-655 nm long-pass 
filter for emission. Time-courses were acquired 10 minutes after 
change in glucose concentration, diazoixide or 1 S-a-glycyrrhetinic 
acid application. Time-courses are displayed as normalized to the 
average fluorescence. 

Data Analysis 

All images were analyzed using custom MATLAB (Mathworks) 
routines or using Slidebook (31). To calculate islet activity, images 
were smoothed using a 5x5 average filter. The variance of pixel 
time-courses was first calculated for a quiescent reference cell; 
manuaUy selected from an area which displayed no fluctuations in 
intensity over time compared with image noise. A pixel was 
considered 'active' if its time-course showed a variance >2 
standard deviations above the variance of the quiescent reference 
cell [1 0,20] . Photobleaching was accounted for through a linear fit, 
and time-courses were rejected if excessive motion artifacts 
occurred. The area of active cells in terms of pixels, was 
determined for each condition and expressed as a fraction of total 
islet pixel area as defined by mean FuraRed fluorescence. GFP 
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regions were defined as liaving a mean pixel fluorescence intensity 
above that measured in GFP~ wild-type cells. The area of active 
cells in GFP"*" regions was expressed as a fi'action of the total GFF*" 
area. Information describing activity is represented as a false-color 
HSV image where Hue is set to 1 (red), Saturation is set to 1 for 
active cells and 0 (no color) for inactive cells, and Value (intensity) 
is set to the average FuraRed fluorescence. Data are presented as 
mean±SEM. For comparison of two means. Student's t-test was 
utilized. For comparison of multiple means, ANOVA was utilized 
along with Tukey's HSD test. 

Boolean Network (Percolation) Model of Suppression 

Bond percolation is a sub-model of percolation theory [31,60] 
which can be used to simulate the islet [14,20,21]. For a lattice of 
nodes (cells) in a given geometry, adjacent nodes are connected 
with a 'coupling probability' p, or not connected with a probability 
{l-p). Connected nodes are considered 'functionally coupled', 
where activity is synchronized at high glucose and suppression 
mediated at low glucose. We implemented simulations of bond 
percolation lattices (figure SI) as previously described [20]. Briefly, 
cubic lattices with alternating node and bond sites (length L= 11) 
were generated. Probabilities were assigned to each bond site with 
a uniform distribution (0 to 1). Neighboring nodes were coupled if 
the bond site probability was less than or equal to the coupling 
probability p. Clusters of coupled nodes were identified and 
potential bond sites removed to establish a matrix of identified 
coupled nodes (figure SI A). 

Coupling-mediated suppression is based on the principle that a 
threshold fraction of non-responsive ('inexcitable') cells can 
suppress all other cells to which they are coupled [17,20]. Within 
a cluster of coupled nodes, if the fraction of inexcitable cells is 
greater than a threshold fraction {Sp), then all cells within the 
cluster are inactive. Experimental studies indicate this threshold is 
<30% [17] and has been modelled to be -15% for MIN6 
aggregates [20] . The probability Pexc defines the fraction of cells 
within the islet that are intrinsically excitable; where in the absence 
of coupUng they would be active. Within each cluster of coupled 
nodes a binomial distribution was used to estimate the probability 
of there being a threshold number of inactive cells to lead to 
suppression. Given a threshold number of inactive cells required 
for suppression (k), a total number of cells in a cluster (m), and the 
fraction of inactive cells (^\ the probability that a coupled cluster is 
active (Pr) is: 



Pr(X = A:) = 



n — k 



where 



{n-k)\k\ 



(2) 



squared minimization. To determine the probability distribution 
for p and Sp, 4000 simulations were run and each simulation was 
separately fitted for p and Sp by chi-squared minimization. 

Coupled Dynamical Model 

The islet model was based on the Cha-Noma fi-cell model 
[29,30], itself based on die Fridlyand P-ceU model [61,62], and 
adapted to include cell-cell coupling. We also included further 
aspects of cell-cell coupling and altered Ka i p channel function. A 
list of parameters used in the model is included in SI Table S 1 . 
The membrane potential (F,) of each P-ceU i is related to the total 
transmembrane current as: 



dVi 

- C — — =IcaV + JtRPM + JsOC + IbNSC + ^KDr + iKCaiSK) , , 

dt (3) 



+ ^KCaiSK) + Ik(A TP) + InuK + InuCu + IpMCA 



Where the kinetics of each current is described in [29,30]. 
To simulate gap junction coupling and a multicellular islet, 
multiple 'cells' were simulated with a coupling current between 

each neighboring cell. The membrane potential for each cell i was 
modified to account for coupling toj neighboring cells: 



-c 



dt 



(4) 



Heterogeneity in coupling was introduced by randomly 
assigning the gap junction conductance g"*';'!, bi^tween cells i and 
J, according to an experimentally measured distribution [unpub- 
lished data], with SD/mean = 70%. To more accurately model |3- 
ceU coupling architecture, random cell lattices were created using 
a position- and availability-based sphere-packing algorithm 
(mean,SD number of cell-cell connections = 5.3,1.7) [63] (figure 
S3). 

Ik(atp) was described in [29,30] as: 



Ik(A TP) = gK(A TP) 'POK(A TPyiV —Vk) 



(5) 



Where gK(ATP) is the open channel conductance and POk(atp) 
represents the mean open probability which is given by: 



0.08 1 



POK(ATP) - 



2[ADP\\ 



([ADP\ 



0.01 



(6) 



0.026 



0.05 J 



The cumulative distribution for k to n, P(X£A), represents the 
probabifity of suppression (sufficient inexcitable cells) in a coupled 
cluster. To obtain the resultant % activity, P(X<A) was averaged 
over all clusters within the islet, weighted by the number of cells n 
in each cluster. '^' normalized to 'n' gives the fraction of 
inexcitable cells required for suppression [Sp). '\-q' gives the 
fraction of excitable cells (Pexc)- 500 simulations were run for each 
P( x( = [0 1] and /) = [0 1], at 0.01 increments for given values of 
Sp. The p and Sp parameters that generated the best fit for the 
simulation mean to experimental data were determined by a chi- 



Endogenous heterogeneity was modelled by randomizing all 
parameters indicated in TableSl between cells about a mean value 
according to a Gaussian distribution with SD/mean = 10%. To 

generate heterogeneity in electrical responses equivalent to 
experimental measurements, the open channel conductance 
gK(ATP) was randomized between cells about a mean value 
according to a Gaussian distribution with SD/mean = 25%. This 
heterogeneity achieves variability in activity that matches exper- 
imental measurements in islets lacking Cx36 [10,14]. 

To model Kir6.2'^'^'^^'^'^'^^^ expression, the open probability 
P0K{ATP) was modified in a proportion (Pexc) of simulated cells: 
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POK(ATP)Mut = yipOKiATP)) + (1 - y){pOMut) (7) 

where y is a constant representing the fraction of ATP-insensitive 
current (pouut = 1)> and was set to 0.5. 

To model diazoxide treatment, the fraction of ATP-insensitive 
current was increased in all cells uniformly according to: 

POKiATP)' = «(P0K(ATP)) + (1 - Cl)(pOMut) (8) 

such that a = 1 represents an untreated islet, and a = 0.5 is 
equivalent to 100% expression of Kir6.2'^^''"''''"-'^^. 

All simulations were initially written and verified in MATLAB, 
then rewritten in C-H- and simulated on the University of 
Colorado JANUS supercomputer. The model was solved using a 
constant time-step Euler integration scheme (Boost C-H- Libraries) 
with 100 [is step-time and 100 ms sampling-time. Rendering 
simulations was performed with Mathematica 9.0 (Wolfram 
Research). 

Supporting Information 

Figure SI Description for how the Boolean network model is 

constructed and how the network parameters P,.^, affects 
network activity. A) Generation of coupled network and 
identification of connected clusters within the network. A 3 
dimensional array of nodes (cells) are generated (in this case an 
L = 11 lattice). Bond sites between neighboring nodes are 
populated according to the coupling probability p, and nodes 
that belong in connected clusters identified. The number of nodes 
in each connected cluster of the network is then recorded. B) 
Histogram for the density of nodes as a function of identified 
cluster sizes for different values oip. Below ~0.25 only very small 
clusters are present. Above ~0.25 a single large connected cluster 
emerges to which the majority of cells belong. Data are averaged 
over 5000 simulations. C) Dependence of cluster activity on the 
cluster size for two different values of Pexo where Sp — 0.15. When 
Pexc<(l-Sp) increasing cluster size leads to reduced activity, 
whereas when Pexc>(l-Sp) increasing cluster size leads to 
increased activity. This explains the respective low and high 
activity in figure IB. Data are averaged over 5000 simulations. 
Overall, the txjupling probability p (representing gap junction 
strength) determines the distribution of connected cluster sizes 
within the islet network. The P^xc value (representing cellular 
excitability) then determines how active each connected cluster is 
depending on its size. The combination of these 2 factors then 
determines the overall network activity. 
(TIF) 

Figure S2 The effect of functional coupUng parameter p on the 
activity of simulated Boolean networks and comparison with 
percolation theory, see [Hraha, T.H., et al., Biophys J, 2014. 
106(1): p. 299-309]. A.) Percent activity of L = 11 lattices was 
simulated for Pexc = 70% for ^ = 0.1 (/)<^c = 0.2488 for cubic 
lattice, i.e. sub-critical coupling) and p = 0.3 (p>Pa i-e. supra- 
critical coupling). B.) Three-dimensional false color maps of 
network activity in representative L = 3 1 lattices for ^ = 0.1 and 
p — 0.3, for Pj.x,; indicated in A. C) Mean size of largest connected 
clusters as a fraction of 3D network in simulated networks as a 
function of ^. D.) Three-dimensional false color maps showing 
distinct 'clusters' of coupled cells in simulated Boolean networks 
with p = 0.2 (p<p,) and^ = 0.3 {p>p,), as indicated in C. Note for 
p = Q.\, the majority of clusters approach single nodes (cells). 
(TIF) 



Figure S3 Penetrance of Kir6.2[^''""'^"'-'QJ and Kir6.2['^]. 
Representative GFP images showing the presence of mutant 
Kir6.2 expressing cells as indicated by GFP co-expression (for 
Kir6.2[^^''"'^"''«J) or a GFP tag (for ¥lr&.2^^\ GFP tagged 
cells are constitutively inactive (inexcitable) in the case of 
Kir6.2'-'^^"^'^'*^"''^'^ and constitutively active for Kir6.2'^'^^^. As 
such Pexc which represents the number of excitable cells is equal to 
100%-%GFP upon Kir6.2f'^^°'''"'^^ expression and is equal to 
%GFP for Kir6.2[^] expression. For Kir6.2[^'''"'''"'''Q^ ex- 
pressing islets, controlled expression is induced through variable 
doses of tamoxifen injections, for which representative images are 
shown for wild-type (GFP = 0%, Pexc =100%), pre-critical 
(GFP =10%, Pexc = 90%), critical (GFP = 20%, Pexc = 80%), and 
post-critical (GFP = 55%, Pexc = 45%), conditions of islet activity. 
Conversely, Kir6.2^'^^^ expressing islets show on average ~70% 
penetrance, such that Pexc = 70%. 
(TIF) 

Figure S4 [Ca^^; in Kir6.2'^'^^^''''^'^^^ expressing and non- 
expressing ceils within the islet as indicated by GFP coexpressioii. 

A) Percent cells showing [Ca'^^]; elevations in GFP positive or GFP 
negative cells as a function of Pexc (100-%GFP). Green diamonds 
indicate those cells expressing GFP whereas black squares indicate 
those cells lacking GFP. B) Mean(±s.e.m.) percent cells showing 
[Ca^""^; elevations in GFP positive or GFP negative cells for data 
binned to wild-type, pre- and post-critical ranges, as determined 
by %GFP. *indicates significant difference (p<0.05) between data 
as indicated. 

(TIF) 

Figure S5 Effects of islet architecture on phase transition 
behavior. A) Schematic representations of the cubic lattice, and 

B) representative example of the quasi-spherical sphere packing 
architectures used for simulating the dynamical oscillator model. 

C) Comparison of simulated islet activity as a function of percent 
excitable cells (Pexc) for cubic and sphere packing architectures in 
the dynamical model for physiological wild-type gap junction 
conductance. 

(TIF) 

Figure S6 Effects of endogenous cellular heterogeneity and 
coupling heterogeneity on phase transition behavior. A) Simulated 

islet activity as a function of percent excitabk- cells (P,.xc) in the 
presence and absentx of heterogeneous distributions of gap 
junction coupling conductance in the dynamical model. B.) 
Simulated islet activity as a function of percent excitable cells 
(Pexc) in the presence and absence of heterogeneous distributions of 
cell physiology parameters in the dynamical model. 
(TIIO 

Figure S7 Dependence of phase transitions in |3-ceU network 
activity on the origin of endogenous P-cell heterogeneity in simulated 
islets. A) Left: Percent cells showing [Ca^^; elevations in simulated 
islet as a function of a uniform increase in the fraction of ATP- 
insensitivity of K^^tp channel activation (a) across cells of the islet. 
Heterogeneity is present in all parameters described in Table S 1 and 
as used elsewhere in this study. Mean simulation data is presented for 
zero gap junction conductance (0 pS) and wild-type gap junction 
conductance (120 pS). Right: activit)' of fully-coupled islet system as a 
function of activity in the uncoupled islet systems which represents the 
excitability of the constituent cells, with heterogeneity present in all 
parameters. For islets lacking gap junction coupling, with zero gap 
junction conductance, the result is trivially linear (dashed). B) As in A 
for heterogeneity solely in Kyip channel activity. C) As in A for 
heterogeneity solely in glycolytic flux. 
(TIF) 
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Figure S8 Mean-field theory analogy of P-cell network activity. 
The excitabiUty of constituent units (i.e. their glucose sensitivities) 
and the resulting network activity takes into account coupling and 
different experimental perturbations. A) Cells of wild-type islets are 
inexcitable at low glucose (2 mM) and all excitable at high glucose 
(20 mM), therefore wild-type islets are res|)ectively fully inactive and 
fully active. Cells expressing Kir6.2'^ N30,ki85Q] glucose- 
insensitive and constitutively inexcitable. When Kir6.2'^'^^°'^'®^^ 
penetrance is <15% (P,.x, >0.85) there are insufficient inexcitable 
cells to suppress global activity, and so coupling leads to inexcitable 
cells being recruited to be active. However, when Kir6.2'^'^^'^'^'^^^ 
penetrance rises above ~15% (Pexc<0.85), global quiescence ensues 
where the majority of normally excitable cells are rendered inactive. 
B) For low diazoxide {<50 |iM), there are fewer excitable cells 
compared to untreated, however if coupling exists then all cells are 
recruited to be artive. However, in the absence of coupling the 
resulting activity is the same as the composition, which results in 
lower activity. In the case of higher diazoxide treatments (> 1 00 [lM) 
the proportion of inexcitable cells exceeds the threshold for 
suppression. If coupling exists then all cells are rendered inactive. 
However in the absence of coupling the resulting activit}' is the same 
as the composition and some cells remain active. C) This mechanism 
can also explain how the islet maintains a robust well-defined 
glucose-stimulated response, but not in the absence of gap junction 
coupling. As glucose is increased more cells become excitable, but in 
the presence of coupling if the excitable fraction is less than a critical 
threshold (e.g. at ~5 mM glucose), all cells are rendered inactive and 
the islet is quiescent. At a glucose level where the excitable fraction 
exceeds the critical threshold (e.g. ~11 mM glucose) all cells are 
recruited to be active. However in the absence of coupling the 
activit)' again reflects the excitability composition within the islet: 
under conditions where normal coupling leads to quiescence greater 
activity occurs, but under conditions where normal coupling leads to 
global activity reduced activity occurs. 
(TIF) 

Movie SI Representative movie of experimentally measured 
[Ca for control conditions. Data is equivalent to that in 
figure 2B, panel I. Note decreases in FuraRed fluorescence are 
equivalent to [Ca^^; increases. Movie played at x20 speed. 
(AVI) 

Movie S2 Representative movie of experimentally measured 
[Ca^-^i for 'pre-criticaf conditions. Data is equivalent to that in 
figure 2B, panel II. Note decreases in FuraRed fluorescence are 
equivalent to [Ca^^j increases. Movie played at x20 speed. 
(AVI) 
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Movie S6 Representative movie of simulated [Ca^"*]; for 'pre- 
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